In [1]:
from netCDF4 import Dataset, num2date
import glob
from mpl_toolkits.basemap import Basemap
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import calendar
from scipy.stats import describe
from scipy import interpolate
from tqdm import tqdm_notebook as tqdm
import oceansdb
In [2]:
def load_netcdf(path = "Argo_South_60"):
files = glob.glob("data/" + path + "/**/*.nc", recursive=True)
lats = []
lons = []
datetimes = []
sst = []
max_depth = []
for f in tqdm(files):
d = Dataset(f)
lat = d.variables["LATITUDE"][:]
mask = lat < -60
lon = d.variables["LONGITUDE"][:]
lats.extend(lat[mask])
lons.extend(lon[mask])
juld = d.variables["JULD"][:]
units = d.variables["JULD"].getncattr('units')
dates = num2date(juld, units, "standard")
datetimes.extend(dates)
try:
sst.extend(d.variables["TEMP_ADJUSTED"][mask,0])
except:
sst.extend(np.full(len(mask), np.nan))
max_depth.extend(np.nanmax(d.variables["PRES_ADJUSTED"], axis=1))
lats = np.array(lats)
lons = np.array(lons)
datetimes = np.array(datetimes)
sst = np.array(sst)
np.save(path + "_lats", lats)
np.save(path + "_lons", lons)
np.save(path + "_dts", datetimes)
np.save(path + "_sst", sst)
np.save(path + "_depth", max_depth)
return lats, lons, datetimes, sst, max_depth
def plot(lats, lons, z = [], title = "Argo profiles south of 60S", cbtitle = "Number of points in bin", vmax=None):
# setup north polar stereographic basemap.
# The longitude lon_0 is at 6-o'clock, and the
# latitude circle boundinglat is tangent to the edge
# of the map at lon_0. Default value of lat_ts
# (latitude of true scale) is pole.
fig = plt.figure(figsize=(15,15))
m = Basemap(projection='spstere',boundinglat=-55,lon_0=180,resolution='l')
m.drawcoastlines()
m.fillcontinents(color='black',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-60, 0, 20))
m.drawmeridians(np.arange(-180, 181, 20), labels=[1,1,0,1])
#m.drawmapboundary(fill_color='aqua')
plt.title("{} {}".format(len(lats), title))
x, y = m(lons, lats)
if len(z) == 0:
hh, locx, locy = np.histogram2d(x, y, bins=100)
# Sort the points by density, so that the densest points are plotted last
z = np.array([hh[np.argmax(a<=locx[1:]),np.argmax(b<=locy[1:])] for a,b in zip(x,y)])
#print(describe(z))
idx = z.argsort()
x, y, z = x[idx], y[idx], z[idx]
#m.imshow(heatmap, interpolation='bicubic', cmap="jet")
m.scatter(x, y, c=z, s=5, alpha=1, cmap="jet", vmax=vmax)
cb = plt.colorbar()
cb.set_label(cbtitle, rotation=270)
plt.show()
def plot_time(dts, title, label):
fig, ax = plt.subplots(1, 1, figsize=(20,10))
for i, dt in enumerate(dts):
k = label[i]
df = pd.DataFrame({k: dt})
df = df.groupby(by=df[k].dt.date).count()
df.plot(ax=ax, style='.', markersize=3)
ax.xaxis.set_major_locator(mdates.YearLocator())
ax.set_title(title)
ax.set_xlabel("Time",fontsize=12)
ax.set_ylabel("Number of profiles",fontsize=12)
#ax.set_xlabel("")
ax.legend()
plt.show()
In [3]:
argo_lats, argo_lons, argo_dts, argo_sst, argo_depth = load_netcdf()
plot(argo_lats, argo_lons, vmax=100)
plot(argo_lats, argo_lons, argo_sst[:len(argo_lats)], title = "Argo SST", cbtitle = "Temperature °C")
plot(argo_lats, argo_lons, argo_depth[:len(argo_lats)], title = "Argo Depth", cbtitle = "Depth (decibar pressure)")
In [4]:
print(min(argo_dts), max(argo_dts))
plot_time([argo_dts], "Argo float dailystamps", ["Argo"])
In [5]:
seal_lats, seal_lons, seal_dts, seal_sst, seal_depth = load_netcdf("seal")
plot(seal_lats, seal_lons, title = "Seal data south of 60S", vmax=100)
plot(seal_lats, seal_lons, seal_sst[:len(seal_lats)], title = "Seal SST", cbtitle = "Temperature °C")
plot(seal_lats, seal_lons, seal_depth[:len(seal_lats)], title = "Seal Depth", cbtitle = "Depth (decibar pressure)")
In [7]:
print(min(seal_dts), max(seal_dts))
plot_time([seal_dts], "Seal data dailystamps", ["Elephant seal"])
In [41]:
all_lats = np.concatenate((seal_lats, argo_lats))
all_lons = np.concatenate((seal_lons, argo_lons))
plot(all_lats, all_lons, title="Argo + seal data south of 60S", vmax=100)
In [8]:
plot_time([seal_dts, argo_dts], "Argo + seal data dailystamps", ["Seal", "Argo"])
In [9]:
def load_netcdf_grid(path = "Argo_South_60", resolution = 1, target_var = "temp", with_depth = False, filter_month = None):
if with_depth:
grid = np.full(shape=[15 * resolution, 360 * resolution, 200], fill_value=np.nan)
else:
grid = np.full(shape=[15 * resolution, 360 * resolution], fill_value=np.nan)
files = glob.glob("data/{}/**/*.nc".format(path), recursive=True)
for f in tqdm(files[:]):
d = Dataset(f)
lat = d.variables["LATITUDE"][:]
mask = (lat < -60) & (lat > -74.5)
if filter_month:
juld = d.variables["JULD"][:]
units = d.variables["JULD"].getncattr('units')
dates = num2date(juld, units, "standard")
datemask = np.array([d.month == filter_month for d in dates])
mask &= datemask
if any(mask):
lat = np.round((np.abs(lat[mask]) - 60) * resolution).astype(int)
lon = d.variables["LONGITUDE"][mask]
lon = np.round((lon + 180) * resolution).astype(int)
lon[lon == 360 * resolution] = 0
if with_depth:
pres = np.round(d.variables["PRES_ADJUSTED"][mask] / 10).astype(int)
temp = d.variables["TEMP_ADJUSTED"][mask]
if target_var == "sal":
try:
sal = d.variables["PSAL_ADJUSTED"][mask]
except:
print("No salinity for {}".format(f))
continue
else:
pres = d.variables["PRES_ADJUSTED"][mask]
temp = d.variables["TEMP_ADJUSTED"][mask, 0]
for x in np.unique(lon):
for y in np.unique(lat):
if with_depth:
ptmask = (lon == x) & (lat == y)
for z in np.unique(pres[ptmask]):
if z>= 200 or np.ma.is_masked(z):
continue
depthmask = pres[ptmask] == z
if target_var == "temp":
mean_for_pt = np.ma.mean(temp[ptmask][depthmask])
elif target_var == "sal":
mean_for_pt = np.ma.mean(sal[ptmask][depthmask])
if not np.isnan(mean_for_pt) and not np.ma.is_masked(mean_for_pt):
grid[y, x, z] = np.nanmean((grid[y, x, z], mean_for_pt))
else:
ptmask = (lon == x) & (lat == y)
if target_var == "n_point":
v = np.sum(ptmask)
if not np.ma.is_masked(v):
grid[y, x] = np.nansum((grid[y, x], v))
elif target_var == "temp":
temps_at_pt = temp[ptmask]
v = np.ma.mean(temps_at_pt)
if not np.ma.is_masked(v):
grid[y, x] = np.nanmean((grid[y, x], v))
elif target_var == "depth":
pres_at_pt = pres[ptmask]
if len(pres_at_pt):
v = np.ma.max(pres[ptmask])
if not np.ma.is_masked(v):
grid[y, x] = np.nanmax((grid[y, x], v))
return grid
def plot_grid(grid, resolution = 1, title="Gridded SST mean for argo data", cbtitle="Temperature °C", vmin=None, vmax=None):
fig = plt.figure(figsize=(15,15))
m = Basemap(projection='spstere',boundinglat=-55,lon_0=180,resolution='l')
m.drawcoastlines()
m.fillcontinents(color='black',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-60, 0, 20))
m.drawmeridians(np.arange(-180, 181, 20), labels=[1,1,0,1])
#m.drawmapboundary(fill_color='aqua')
plt.title("{} at {}° resolution".format(title, 1/resolution))
x = np.arange(-180, 180.1, 1/resolution)
y = np.arange(-60, -75.1, 1/-resolution)
x, y = np.meshgrid(x, y)
x, y = m(x, y)
m.pcolormesh(x, y, grid, cmap="jet", vmin=vmin, vmax=vmax)
cb = plt.colorbar()
cb.set_label(cbtitle, rotation=270)
plt.show()
def interp_nans(grid, method='linear'):
x = np.arange(0, grid.shape[1])
y = np.arange(0, grid.shape[0])
mask = np.isnan(grid)
xx, yy = np.meshgrid(x, y)
x1 = xx[~mask]
y1 = yy[~mask]
newarr = grid[~mask]
interp = interpolate.griddata((x1, y1), newarr.ravel(), (xx, yy), method=method)
if method == "linear":
interp = interp_nans(interp, "nearest")
return interp
In [46]:
resolution = 2
argo_pt_grid = load_netcdf_grid(target_var="n_point", with_depth = False, resolution = resolution)
seal_pt_grid = load_netcdf_grid(path="seal", target_var="n_point", with_depth = False, resolution = resolution)
In [47]:
combined_pt_grid = np.nansum((argo_pt_grid, seal_pt_grid), axis=0)
combined_pt_grid[np.isnan(argo_pt_grid) & np.isnan(seal_pt_grid)] = np.nan
In [48]:
plot_grid(argo_pt_grid, title="Number of points in gridded argo dataset", cbtitle="Number of points", vmax=50, resolution = resolution)
plot_grid(seal_pt_grid, title="Number of points in gridded seal dataset", cbtitle="Number of points", vmax=50, resolution = resolution)
plot_grid(combined_pt_grid, title="Number of points in gridded argo+seal datasets", cbtitle="Number of points", vmax=100, resolution = resolution)
In [222]:
resolution = 4
In [223]:
argo_temp_grid = load_netcdf_grid(resolution = resolution)
seal_temp_grid = load_netcdf_grid("seal", resolution = resolution)
In [224]:
combined_temp_grid = np.nanmean((argo_temp_grid, seal_temp_grid), axis=0)
In [225]:
plot_grid(argo_temp_grid, resolution = resolution)
plot_grid(seal_temp_grid, title = "Gridded SST mean for seal data", resolution = resolution)
plot_grid(combined_temp_grid, title = "Gridded SST mean for argo+seal data", resolution = resolution)
In [226]:
interp = interp_nans(combined_temp_grid)
In [227]:
plot_grid(interp, resolution = resolution, title="Interpolated gridded SST mean for argo + seal data")
In [155]:
argo_max_depth = load_netcdf_grid(target_var = "depth", resolution = resolution)
seal_max_depth = load_netcdf_grid("seal", target_var = "depth", resolution = resolution)
In [159]:
plot_grid(argo_max_depth, resolution = resolution, title="Argo gridded max depth", cbtitle="Depth (dbar)")
plot_grid(seal_max_depth, resolution = resolution, title="Seal gridded max depth", cbtitle="Depth (dbar)", vmax=2000)
In [281]:
for month in range(1, 13):
argo_temp = load_netcdf_grid(filter_month = month, with_depth=True)
seal_temp = load_netcdf_grid("seal", filter_month = month, with_depth=True)
combined_temp = np.nanmean((argo_temp, seal_temp), axis=0)
np.save("argo_temp_grid_{}".format(month), argo_temp)
np.save("seal_temp_grid_{}".format(month), seal_temp)
np.save("combined_temp_grid_{}".format(month), combined_temp)
interp_temp = interp_nans(combined_temp[:,:,0])
argomean = np.round(np.nanmean(argo_temp), 2)
sealmean = np.round(np.nanmean(seal_temp), 2)
combmean = np.round(np.nanmean(combined_temp), 2)
interpmean = np.round(np.nanmean(interp_temp), 2)
title = "Argo 0-10 decibar mean for {}. Mean across all points = {}".format(calendar.month_abbr[month], argomean)
plot_grid(argo_temp[:,:,0], title=title)
title = "Seal 0-10 decibar mean for {}. Mean across all points = {}".format(calendar.month_abbr[month], sealmean)
plot_grid(seal_temp[:,:,0], title=title)
title = "Argo+Seal 0-10 decibar mean for {}. Mean across all points = {}".format(calendar.month_abbr[month], combmean)
plot_grid(combined_temp[:,:,0], title=title)
title = "Interpolated Argo+Seal 0-10 decibar mean for {}. Mean across all points = {}".format(calendar.month_abbr[month], interpmean)
plot_grid(interp_temp, title=title)
In [ ]:
resolution = 2
try:
argo_temp_grid_withdepth = np.load("argo_temp_grid_withdepth.npy")
seal_temp_grid_withdepth = np.load("seal_temp_grid_withdepth.npy")
combined_temp_grid_withdepth = np.load("combined_temp_grid_withdepth.npy")
except:
argo_temp_grid_withdepth = load_netcdf_grid(with_depth=True, resolution = resolution)
seal_temp_grid_withdepth = load_netcdf_grid("seal", with_depth=True, resolution = resolution)
combined_temp_grid_withdepth = np.nanmean((argo_temp_grid_withdepth, seal_temp_grid_withdepth), axis=0)
np.save("argo_temp_grid_withdepth", argo_temp_grid_withdepth)
np.save("seal_temp_grid_withdepth", seal_temp_grid_withdepth)
np.save("combined_temp_grid_withdepth", combined_temp_grid_withdepth)
In [91]:
dbars = [30, 50, 100, 150, 199]
vmax = 5
vmin = -3
for db in dbars:
dbs = "from {}-{} dbar".format(db*10, (db+1)*10)
plot_grid(argo_temp_grid_withdepth[:,:,db], title="Gridded argo temperature data " + dbs, resolution=resolution, vmin=vmin, vmax=vmax)
plot_grid(seal_temp_grid_withdepth[:,:,db], title="Gridded seal temperature data " + dbs, resolution=resolution, vmin=vmin, vmax=vmax)
plot_grid(combined_temp_grid_withdepth[:,:,db], title="Gridded seal+argo temperature data " + dbs, resolution=resolution, vmin=vmin, vmax=vmax)
interp = interp_nans(combined_temp_grid_withdepth[:,:,db])
plot_grid(interp, title="Interpolated gridded seal+argo temperature data " + dbs, resolution=resolution, vmin=vmin, vmax=vmax)
In [7]:
fig = plt.figure(figsize=(15,15))
m = Basemap(projection='spstere',boundinglat=-55,lon_0=180,resolution='l')
#m.drawcoastlines()
#m.fillcontinents(color='black',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-60, 0, 20))
m.drawmeridians(np.arange(-180, 181, 20), labels=[1,1,0,1])
#m.drawmapboundary(fill_color='aqua')
m.etopo()
plt.title("Bathymetry south of -60S")
Out[7]:
In [220]:
def plot_cross_slope(grid, title="Temperature at 179° W", lon=-179, cbtitle="Temperature °C", resolution = 1):
fig, ax = plt.subplots(1, 1, figsize=(20,10))
y = np.arange(-60, -75, -1/resolution)
z = np.arange(0, 2000, 10)
levels = np.arange(np.nanmin(grid), np.nanmax(grid), .01)
im = ax.contourf(y, z, grid.T, cmap="jet", levels = levels, extend="both")
ax.set_xlabel("° Latitude")
ax.set_ylabel("Depth (dbar)")
ax.set_title("{} at {} resolution".format(title, 1 / resolution))
cb = fig.colorbar(im)
cb.set_label(cbtitle, rotation=270)
db = oceansdb.ETOPO()
y = np.arange(-60, -75, -.01)
h = -db['topography'].extract(lat=y, lon=lon)["height"]
#ax.plot(y, h)
ax.fill_between(y, 5000, h, color='black')
ax.set_ylim(2000, 0)
ax.set_xlim(-60, -74.5)
plt.show()
filtered = combined_temp_grid_withdepth[:,2,:]
interp = interp_nans(filtered)
plot_cross_slope(filtered, resolution = resolution)
plot_cross_slope(interp, title="Interpolated Temperature at 179° W", resolution = resolution)
In [94]:
try:
argo_sal_grid_withdepth = np.load("argo_sal_grid_withdepth.npy")
seal_sal_grid_withdepth = np.load("seal_sal_grid_withdepth.npy")
combined_sal_grid_withdepth = np.load("combined_sal_grid_withdepth.npy")
except:
argo_sal_grid_withdepth = load_netcdf_grid(with_depth=True, target_var="sal", resolution = resolution)
seal_sal_grid_withdepth = load_netcdf_grid("seal", with_depth=True, target_var="sal", resolution = resolution)
combined_sal_grid_withdepth = np.nanmean((argo_sal_grid_withdepth, seal_sal_grid_withdepth), axis=0)
np.save("argo_sal_grid_withdepth", argo_sal_grid_withdepth)
np.save("seal_sal_grid_withdepth", seal_sal_grid_withdepth)
np.save("combined_sal_grid_withdepth", combined_sal_grid_withdepth)
In [143]:
dbars = [30, 50, 100, 150, 199]
vmax = 35
vmin = 33.5
for db in dbars:
dbs = "from {}-{} dbar".format(db*10, (db+1)*10)
interp = interp_nans(combined_sal_grid_withdepth[:,:,db])
plot_grid(interp, title="Interpolated gridded seal+argo salinity data " + dbs, cbtitle="Salinity", resolution=resolution, vmin=vmin, vmax=vmax)
In [153]:
def dens_smow(T):
'''
Function calculates density of standard mean ocean water (pure water) using EOS 1980.
INPUT: T = temperature [degree C (ITS-90)]
OUTPUT: dens = density [kg/m^3]
'''
a0 = 999.842594;a1 = 6.793952e-2;a2 = -9.095290e-3;a3 = 1.001685e-4;a4 = -1.120083e-6;a5 = 6.536332e-9
T68 = T * 1.00024
dens = (a0 + (a1 + (a2 + (a3 + (a4 + a5*T68)*T68)*T68)*T68)*T68)
return dens
def dens0(S, T):
'''
Function calculates density of sea water at atmospheric pressure
USAGE: dens0 = dens0(S,T)
DESCRIPTION:
Density of Sea Water at atmospheric pressure using
UNESCO 1983 (EOS 1980) polynomial.
INPUT: (all must have same dimensions)
S = salinity [psu (PSS-78)]
T = temperature [degree C (ITS-90)]
OUTPUT:
dens0 = density [kg/m^3] of salt water with properties S,T,
P=0 (0 db gauge pressure)
'''
assert S.shape == T.shape
T68 = T * 1.00024
# UNESCO 1983 eqn(13) p17.
b0 = 8.24493e-1;b1 = -4.0899e-3;b2 = 7.6438e-5;b3 = -8.2467e-7;b4 = 5.3875e-9
c0 = -5.72466e-3;c1 = +1.0227e-4;c2 = -1.6546e-6
d0 = 4.8314e-4
dens = (dens_smow(T) + (b0+ (b1+(b2+(b3+b4*T68)*T68)*T68)*T68)*S + (c0+(c1+c2*T68)*T68)*S*np.sqrt(S)+ d0*(S**2) )
return dens
In [166]:
dens = dens0(combined_sal_grid_withdepth, combined_temp_grid_withdepth)
dbars = [30, 50, 100, 150, 199]
vmax = 1028
vmin = 1026.5
for db in dbars:
dbs = "from {}-{} dbar".format(db*10, (db+1)*10)
interp = interp_nans(dens[:,:,db])
print(np.nanmin(interp), np.nanmax(interp))
plot_grid(interp, title="Interpolated gridded seal+argo density data " + dbs, cbtitle="Density", resolution=resolution, vmin=vmin, vmax=vmax)
In [219]:
plot_cross_slope(dens[:,2,:], title="Density at 179° W", cbtitle="Density", resolution=resolution)
plot_cross_slope(interp_nans(dens[:,2,:]), title="Interpolated density at 179° W", cbtitle="Density", resolution=resolution)
In [ ]: